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ABSTRACT 
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i 

''Solutions  of  the  ray  propagation  equation  are  obtained 
for  various  boundary  conditions,  assuming  vertical  incidence 
in  plane,  parallel,  equal-travel-time  layers.  The  solutions 
are  examined  in  noth  the  time  and  frequency  domains  and 
certain  properties  derived.  A  complete  discussion  of 
frequency-domain  synthesis  techniques  is  given  in  connection 
with  a  treatment  of  the  absorption  problem.  FORTRAN  programs 
are  given  which  compute  any  of  the  solutions  in  either  the 
frequency  or  time  domains,  with  or  without  absorption.  The. 
theory  and  programs  are  applied  to  the  problem  of  source  depth 
determination,  and  it  is  shown  that  the  method  of  pP  spectral 
nulls  is  somewhat  unreliable. 
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INTRODUCTION 


Considerable  work  on  ray  propagation  in  multilayered  media 
has  been  reported  in  the  literature  (see  Claerbout,  1968,  for 
references).  Most  of  it  has  been  concerned  with  calculating 
near-source  seismograms  for  vertically  incident  waves  in  loss¬ 
less  media.  Claerbout  summarized  previous  results  for  this  case 
in  the  concise  notation  of  Sherwood  and  Trorey  (1965),  and 
included  in  his  paper  a  computer  program  to  synthesize  reflection 
seismograms  in  the  time  domain.  He  also  derived  the  solution  for 
the  transmission  seismogram  due  to  a  source  buried  in  the  lower 
half-space.  Landers  and  Claerbout  (1969)  applied  the  half-space 
transmission  solution  to  a  study  of  crust/upper  mantle  models 
proposed  by  Aki.  Frasier  (1970)  has  developed  solutions  for  non¬ 
vertical  ly-i  nci  dent  P  and  SV  waves  in  non-absorpti ve  media  which 
are  closely  analogous  to  Claerbout's. 

A  totally  satisfying  treatment  of  the  absorptive  case  has 
not  been  given,  to  the  author's  knowledge.  Trorey  (1962)  devised 
a  time-domain  solution  using  a  non-real i zable  linear  absorption 
law.  For  practical  reasons  he  was  forced  to  lump  his  absorption 
into  a  small  number  of  constant-Q  bands.  His  procedure  seemed 
cumbersome  and  no  attempt  was  made  by  this  author  to  duplicate 
it.  Instead,  a  frequency-domain  approach  has  been  utilized. 
Several  authors  have  expressed  the  fear  that  this  wou)d  lead 
to  large  aliasing  errors,  but  such  has  not  proved  to  be  the 
case.  A  discussion  of  practical  difficulties  is  given.  Besides 
its  simplicity,  the  frequency-domain  calculation  has  the 
additional  advantage  that  there  is  no  particular  problem  in 
making  one's  absorption  law  realizable.  Sherwood  and  Trorey 
(1965)  gave  essentially  a  physical  argument  that  the  delayed 
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transmission  seismogram  must  be  minimum  phase  for  a  minimum- 
phase  absorption  law.  A  mathematical  proof  is  presented  here. 

A  discussion  of  dispersion  is  also  given. 

Since  this  work  was  completed,  a  paper  has  appeared 
(Jensen  and  Ellis,  1970)  in  which  the  authors  obtain  solutions 
for  non-vertical  incidence  plus  absorption,  using  linear  system 
theory.  They  calculate  spectra,  not  seismograms,  and  they  do 
not  include  dispersion.  No  programs  are  given.  An  alternative 
approach  to  the  general  problem  would  be  to  apply  the  techniques 
developed  here  to  the  solutions  given  by  Frasier. 

The  emphasis  in  the  present  work  has  been  placed  on  obtaining 
algorithms  which  are  fast,  accurate,  and  concise,  and  therefore 
useful  for  large-scale  model  studies.  The  first  part  of  the  paper 
is  devoted  to  obtaining  lossless  transmission  solutions  for 
various  boundary  conditions  of  interest  (assuming  vertical 
incidence  in  plane,  parallel,  equal-travel-time  layers).  A 
computer  program  is  included  to  do  these  calculations  (exactly) 
in  the  time  domain.  The  second  part  considers  the  practical 
aspects  of  frequency-domain  absorption  calculations.  FORTRAN 
programs  are  included  to  do  all  the  cases  treated  in  the  first 
part  in  the  frequency  domain,  including  frequency-  and  depth- 
dependent  absorption.  Versions  both  with  and  without  dispersion 
are  given. 

One  application  of  the  theory  and  programs,  to  source  depth 
determination,  is  discussed  in  detail.  Sources  of  error  in 
determinations  by  P-wave  spectral  nulls  or  cepstral  analysis  are 
investigated  both  theoretically  and  experimentally.  It  is  shown, 
in  particular,  that  inhomogeneities  near  the  source  can  produce 
large  shifts  in  the  null  frequencies  from  the  values  predicted 
using  a  simple  echo  model.  Further,  it  appears  that  there  is  no 
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really  practical  way  to  circumvent  this  effect.  Thus  a  limit 
can  be  placed  on  the  reliability  of  this  method. 

Certain  other  applications  of  the  programs  are  mentioned 
briefly. 


THE  LOSSLESS  CASE 


Propagation  equation 

For  up-  and  downgoing  displacement  waves  U  and  D,  waves 
at  the  surface  are  related  to  those  at  depth  by  the  "propagator 
matri x" 


u 

1 

zkF(l/z) 

zkG(l/z) 

U 

D 

_  _ 

surface  w^t 

G(z) 

F(z) 

D 

Derivation  of  (1)  is  given  in  detail  by  Claerbout  and  will  not 
be  repeated  here.  (Note  that  for  convenience  we  factor  out  a 
transmission  factor  IT  t  s  TTk  =  1t..)  F  and  G  are  polynomials  in 
z,  obtained  by  taking  products  of  k  "layer  matrices" 


U 

=  !_ 

z 

zr 

U 

D 

wt 

j 

r 

1 

D 

j  +  1 


(2) 


where  z  =  w ^  -  e"^1  represents  a  unit  delay  operator,  and  x 
is  the  two-way  travel-time  across  each  layer  (the  same  for  all 
layers)  and  equals  the  sampling  interval.  The  (displacement) 
reflection  coefficients  r  are  defined  at  each  interface  by 


r  =  (pava  -  Pbvb)/(Pava  +  Pbvb) 
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for  densities  p  and  velocities  v,  £  above  and  below;  and 
the  (downgoing)  transmission  coefficients  t  are  related  to 
the  reflection  coefficients  by  t  =  1  +  r.  Equation  (1)  may  be 
conveniently  viewed  in  either  the  time  or  frequency  domains. 

In  this  paper  we  consider  only  problems  in  which  at  least 
one  free  surface  is  present,  and  it  turns  out  that  in  this  case 
the  quantity  of  interest  in  (1)  is  not  F  or  6,  but  the  combi¬ 
nation 

A ( z )  =  F(z)  -  zkG ( 1 / z )  (4) 

(A  is  normalized  and  is  related  to  Claerbout's  M  by  A(z)  = 

n.t.M(z).)  A  recursion  for  A  may  be  developed  by  multiplying 
the  propagator  for  k-1  layers  by  an  additional  layer  matrix 


“  - 

“ 

zk-1 F ( 1 / z )  zk_1 G ( 1 /z ) 

2  zrk 

G(z)  F(z) 

rk  1 

_  — , 

rkzkF(l/z)  +  zk_1  G ( 1 / z ) 
rkzG(z)+F(z) 

from  which 


A(k)(z)  *  [rkzG(z)  +  F ( z ) ]  -  [rkzkF(l/z)  +  zk-,G(l/z)] 


by  definition  (4),  or 


(k) 


(z) 


A(  k"'1  J  ( z )  - 


,zkA^k‘] )(l/z) 


(5) 


Claerbout  gives  essentially  this  expression.  In  the  time-domain 
this  is  equivalent  to  the  recursion 


A 


(k) 

1 


(k-1) 


•  it 


1 


A(K) 

Aj 


rkA 


(k-1) 
k-j  +  2 


j  =  2 ,  ....  k 


A 


(k) 

k+1 


-r 


k 


(obtained  by  putting  A^(z)  =  A-|  +  A2z  +  AgZ2  +  ...  +  Ak+1z 
into  (5)  and  identifying  coefficients  of  powers  of  z). 


Half-space 

For  the  case  of  an  observer  located  at  a  free  surface  and 
an  impulsive  source  buried  in  a  half-space  below  the  layers, 
the  boundary  conditions  are  a  perfect  reflection  of  the  observed 
seismogram  X ( z )  at  the  surface,  a  1  coming  up  from  below,  and 
some  function  P(z)  returned  into  the  half-space.  The  propagation 
equation  (1)  is  thus 


X 

v  1 

zkF(l/z) 

zkG(l/z) 

1 

X 

wk  nt 

G  ( z ) 

F  ( z ) 

P 
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Inverting  and  solving  for  X  yields 


X  ( z )  =  wkrr|[=1(l-r,)/A(z)  (7) 

(The  determinant  in  (6)  is  zkF(l/z)F(z)-zkG(l/z)G(z)  =  zkTT1-  ( 1  -r? ) , 
from  (2).)  This  is  Claerbout's  solution.  The  seismogram  due  to 
an  arbitrary  source  function  S(z)  may  be  obtained  by  convolving 
the  impulse  response  with  S. 

To  actually  compute  (7)  we  first  compute  A  by  means  of  the 

recursion  (5)  and  then  invert.  Note  that  in  the  programs  we  do 

our  convolutions  (multiplications)  and  deconvolutions  (divisions) 

in  the  time  domain  for  accuracy.  In  the  timedomain,  convolution 

2 

and  deconvolution  are  N  processes  (require  a  time  proportional 
2 

to  N  for  N  elements),  but  since  the  layer  recursion  (5)  is  also 

2  k 

N  ,  little  is  lost  by  doing  this.  Note  also  that  the  factor  w 
in  (7)  corresponds  to  the  initial  delay  of  the  first  arrival, 
but  this  is  ignored  in  the  programs  and  output  commences  with 
the  first  point. 

Source  on  far  surface 

If  both  observer  and  source  are  located  at  free  surfaces, 
boundary  conditions  are  as  shown  in  Figure  1. 

SURFACE 

Tx  Jx 


♦  hr _ +R 

SURFACE 

Figure  1.  Source  on  far  surface 
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R  here  is  some  unknown  function.  The  propagation  equation 
for  this  case  is 


X 

1 

zkF( 1/z) 

z  k  6  ( 1  /  z ) 

1  +  R 

X 

w  t 

G(z) 

F(z) 

R 

Solving  for  X, 


X(z )  [  A  ( z )  -  zkA(  1  /z )  ]  =  wk  TT.O-r.)  (8) 

Inspection  of  (5)  and  (7)  shows  that  the  same  answer  could  be 
obtained  from  the  half-space  solution  by  adding  an  extra  layer 
with  r  =  1  (an  obvious  result),  provided  that  the  layer  is 
not  included  in  the  transmission  factor  TT ^ ( 1 -r^ ) .  Writing 
A*(z)  =  A(z)  -  zkA(l/z),  (a  definition  chat  will  be  used  through¬ 
out  this  paper),  (8)  becomes 


X  ( z )  =  wkTT(l-r.)/A*(z) 


(8*) 


The  inverse  wavelet  in  this  case  is  antisymmetric  about  its 
midpoint.  This  implies,  among  other  things,  that  half  the 
information  about  the  layers  is  lost,  that  is,  one  cannot  go 
from  the  seismogram  back  to  the  reflection  coefficients  of 
the  layers  in  this  case.  Another  way  to  see  this  is  as  follows. 
Inverting  z  in  (5),  multiplying  by  rkzk  and  adding  to  A  gives 
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(9) 


( 1  “rk  > 


A ( k" 1 ^  ( z ) 


M 


(z)  + 


rkzk  A(k)(l/z) 


which  is  the  inverse  recursion  that  allows  us  to  extract  the 

reflection  coefficients  from  the  inverse  wavelet.  In  (9)  we 

( k ) 

remove  layers  one  at  a  time,  where  at  each  step  rk  =  -Ak+1  is 
the  last  boundary.  Evidently  (9)  blows  up  when  any  rk  =  +  1 . 

This  result  has  an  interesting  physical  interpretation. 
Consider  the  propagator  of  an  inverted  layer  set.  If  the 
normal  case  is  written 


U' 


D' 


=  Q 


U 

D 


we  write  the  inverted  case 


“  ' 

D 

II 

-Ol 

D' 

U 

U' 

by  inspection  of  Figures  2  and  3. 


Figure  2.  Normal  case 
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Figure  3.  Layers  inverted 


Comparison  with  (1)  gives 


S(z) 


1 

wkn(l-r.) 


zkF(l/z) 

-zkG(l/z) 


-6(z) 

F(z) 


(10) 


This  can  also  be  obtained  from  (1)  by  a  time  reversal,  z  1/z, 
plus  a  matrix  inversion,  i.e.  Q(z)  =  Q  ( 1 / z ) .  Writing 
A(z)  =  p(z)  +  g(z),  by  analogy  with  (4),  we  have 

J*(z)  =  A(z)  -  zk  A(l/z) 

=  [F( z)  +  G( z ) ]  -zk[F(l/z)  +  G ( 1 / z ) ] 

=  A(z)  -  zkA(l/z) 


or 


A  * ( z )  =  A*(z) 


(ID 


Thus  the  antisymmetric  combination  is  invariant  under  an 
inversion  of  the  layers.  Comparing  this  with  (8‘),  we  see 
that  interchanging  source  and  receiver  gives  the  same  seismo 
gram,  except  for  a  scale  factor  (reciprocity),  i.e.  one 
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cannot  tell,  from  the  seismogram,  which  side  of  the  earth  is 
which. 

In  principle,  a  reflection  coefficient  slightly  less  than 

unity  removes  the  difficulty  associated  with  (9),  but  in 

2 

practice  the  requirement  is  that  1-r  must  not  be  so  small  as 
to  make  the  recursion  unstable.  In  any  case,  the  actual  reflec¬ 
tion  coefficient  of  the  distant  free  surface  is  just  the  ratio 
of  the  last  to  first  points  of  the  inverse  wavelet.  The 
presence  of  absorption  will  change  these  results  drastically. 

In  fact  the  inverse  problem  is  then  not  do-able  even  in  the 
half-space  case  (at  least  not  using  the  techniques  discussed 
here),  essentially  because  one  is  required  to  obtain  2k  items 
of  information  from  only  k  items  of  data. 

Contained  source 

Of  some  interest  is  the  case  of  a  source  buried  in  a 
layer  stack  terminated  at  both  ends  by  free  surfaces.  This 
case  may  be  easily  treated  by  coupling  half-space  solutions 
back  to  back  (Figure  4)„ 


SURFACE 


W~  X? 


SURFACE 


Figure  4,  Buried  source 
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We  assume  the  source  to  be  located  in  a  trivial  (r  =  0)  inter¬ 
face  common  to  the  two  layc"  stacks.  The  transmitted  wave  due 
to  a  source  B(z)  is,  by  (7), 


X(z)  =  wk  TT.O-r.)  B(z)/A(z) 


(7*) 


The  wave  returned  from  the  half-space  may  be  obtained  by 
solving  (6)  for  P 


P ( z )  =  w k  A ( 1  / z )  X(z)  /TT|  ( 1  - r i ) 
=  z k  B ( z )  A ( 1 /z )  / A ( z ) 


02a) 

(12b) 


Primes  will  be  used  to  denote  quantities  associated  with  the 
distant  stack.  From  Figure  4,  the  coupling  relations  are 


B ( z )  =  S(z)  -  P'(z) 

(13) 

B'(z)  =  S(z)  -  P(z) 

(13') 

for  a  spatially  symmetric  source  S(z).  The  minus  signs  in  (13) 
are  necessary  to  take  account  of  the  change  in  reference 
direction  of  the  (displacement)  waves  between  the  two  layer 
sets.  Substituting  (12a)  into  (13'),  (13‘)  into  (12b)  (primed), 
and  the  result  into  (13)  gives 
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B ( z )  =  S(z)  -  zk'  [S ( z ) - A(  1  /z ) X ( z )  ]  ililM 

TT  ( 1  —  r  - )  A 1  ( z ) 

Using  (7')  and  collecting  terms 

[A( z ) A 1 ( z )  -  zk+k'  A(l/z)A'(l/z)]  X ( z ) 

(14) 

=  wk7T  (1-r.)  [  A '  (z )  -  zk  A 1  ( 1  /  z )  ]  S  ( z ) 
i  =1  1 

Note  that  the  filter  on  the  LHS  of  (14)  is  antisymmetric  and 
depends  on  all  the  layers,  whereas  the  filter  on  the  RHS, 
which  is  also  antisymmetric,  depends  only  on  the  layers  lying 
between  the  source  and  the  far  surface. 


A  comparison  w'th  the  far-surface  solution  (8)  is  tempting, 
but  in  that  case  one  reference  direction  was  used  for  all 
layers,  whereas  in  this  case  we  have  decomposed  the  stack  into 
two  opposing  substacks.  To  make  the  comparison,  we  should  write 
the  propagator  for  the  whole  stack  in  terms  of  the  two  substacks. 
This  is  just  the  product  of  the  propagator  for  the  near  stack 
and  the  prr  ’ator  for  the  far  stack  inverted,  i.e.  QT  =  Q  0  . 

From  (1 )  and  ,  )  thi s  is 


Q 


T 


1 

wkirt 


zkF(l/z)  zkG(l/z) 
G(z)  F(z) 


v*7f  ( 1  -r  i ) 


zk ' F 1 ( 1 /z ) 
-z k  '  G 1 ( 1 / z ) 


-G'(z) 

F'(z) 


1 _ 

wk+k,lT(l  +  ri  )7T(l-r! ) 


-G'(z)zkF(l/z)+F'(z)zkG(l/z) 
-G* (z)G(z)+F* (z)F(z) 
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from  which 


At(z)  =  F'(z)  A ( z )  +  G' (z)zkA(l/z) 
by  (4),  and 

*  kT 

At  5  At(z )-z  t  Aj( 1/z) 

=  A ( z ) A ' (z)-zk+k'  A(l/z)A'(l/z)  =  ( AA ' )* 


(15) 


Putting  this  result  in  (14)  gives 


A*(z)  X  ( z )  =  wk7Tk’1(l-r.)  A'*  (z)  S(z) 
1  i  =  l  1 


(14*) 


Comparison  of  (14*)  with  ( 8 * )  shows  that  the  only  difference  is 
in  the  RHS  of  (14* ),  that  is,  one  obtains  the  same  seismogram 
from  a  source  S$  located  at  the  far  surface  as  from  a  buried 
source  Sb  provided 


Ss(z) 


TT(l+r') 


Sb(z) 


(15) 


(The  change  in  the  sign  of  the  reflection  coefficients  arises 
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because  we  have  considered  the  primed  stack  to  be  ordered 
from  the  far  surface  down.  Actually,  theorem  (11)  says  that 
we  may  calculate  A’*  from  the  reversed  stack  as  well,  in 
order  to  be  consistent  with  the  calculation  of  A  -p.  This  is 
done  in  the  programs).  Thus  (15)  represents  a  "source-burying 
operator"  which  describes  the  effect  of  burying  the  source 
in  terms  of  the  layers  between  the  far  surface  and  the  source. 

Remembering  that  A’*  is  antisymmetric  about  its  midpoint, 
which  represents  the  depth  of  the  source  in  travel-time 
(sampling-rate)  units,  we  might  hope  to  use  this  fact  in 
source  depth  determination.  If  one  could  isolate  A'  by  some 
means,  then  the  least-squares  point  of  antisymmetry  would 
give  the  depth.  For  example,  if  a  seismogram  Xg  were  available 
from  an  identical  source  lying  along  the  same  path  at  the 
surface,  then  comparison  with  the  buried-source  seismogram 
Xb  would  give 


VXs  ~  A'* 

ignoring  scale  and  delay  factors.  This  method  suffers  from 
certain  obvious  disadvantages.  A  more  direct  approach  would 
be  to  look  for  zeros  of  Xfa.  This  has  been  done  in  the  past 
(cepstrum  analysis)  and  suffers  principally  from  the  effect 
of  the  source  window  S(z).  A  study  of  the  limitations  of  the 
method  is  presented  below  under  Applications.  An  alternative 
derivation  of  the  buried-source  operator  given  in  the  next 
section  shows  that  the  result  is  general. 
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Source-burying  operator: 

An  alternative  derivation  of  the  source-burying  operator 
can  be  given  by  comparing  (distant)  half-space  solutions  for  a 
far-surface  source  and  requiring  that  the  incoming  and  outgoing 
waves  be  identical  for  the  two  cases  (Figures  5  and  6). 


© 

t  H  ffcS, 

te  It 

Figure  5.  Surface  source 


fx  Tx 

u  |  *' 

Figure  6.  Buried  source 


For  a  surface  source,  the  propagation  equation  (1)  is 


R 

1 

zkF(l/z) 

zkG(l/z) 

B 

R  +  S$ 

wknt 

G  ( z ) 

F(z) 

T 
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Solving  for  T  gives 


wkn(l+r.)  zkA( 1 /z ) 

T  =  - —  S  +  -  B 

A ( z )  A(z) 


(16) 


The  solution  for  the  buri ed-source  case  is  given  by  (12b)  with 
source  B-S  (the  minus  being  necessary  to  keep  the  reference 
direction  consistent  with  the  previous  case),  i.e. 


P  =  [zkA(l/z)  /  A(z)]  (B-S) 


Then 


P  +  S 


A(z)-zKA(,l,/z 1 
A  ( z ) 


S  + 


zKA(1/z) 

A(z) 


(17) 


Identifying  T  in  (16)  with  P  +  S  in  (17)  gives 


s  =  AUl-Z^Q/z)  s 
s  wkn(l+r.) 


(15) 


as  before. 

In  the  buried-source  solutions  we  have  so  far  considered 
only  the  case  of  a  source  located  within  a  homogeneous  layer.  If 
the  source  is  located  in  a  reflecting  interface,  we  can  modify 
the  above  derivation  by  explicitly  considering  the  reflections 
at  the  additional  interface.  If  this  boundary  has  reflection 
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coefficient  c,  instead  of  (15)  we  obtain 


S  =  A(z)  -  (l+2c)  zkA( 1 / z )  s 
s  w^d  +  r.) 


(18) 


This  resembles  the  previous  result  with,  however,  a  surface 
reflection  coefficient  r  =  l+2c.  (This  can  easily  be  proved 
by  inverting  the  stacks  in  Figures  5  and  6  and  going  through 
the  calculation  assuming  a  surface  reflection  coefficient  r 
other  than  unity). 

We  can  also  modify  the  previous  derivation  to  include  the 
case  of  a  spatially  asymmetric  source.  If  in  Figure  6  we 
replace  the  upcoming  (toward  the  surface)  part  of  the  source 
function  by  S'(z),  then  it  is  easy  to  show  that  (15)  becomes 


s  .  A(Z)S(Z)  -  ZkA(l/z)S'  (z)  09 

s  wkn(l+r.) 

If  S(z)  and  S'(z)  have  identical  time  behavior  but  differ  by 
a  scale  factor  (for  example,  an  earthquake  with  an  asymmetrical 
radiation  pattern)  then  (19)  will  have  the  form  of  (18).  The 
programs  have  not  been  written  to  handle  these  cases,  but  it 
would  be  easy  to  modify  them  to  do  so,  simply  by  redefining 
the  value  of  the  terminal  reflection  coefficient  between  calls 
to  the  inner  computation  routine  (RCTOA). 

Buried-recei ver  operator 

The  results  we  have  obtained  have  assumed  the  receiver  to 
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be  located  at  the  surface.  A  "buri ed-recei  ver  operator", 
analogous  to  the  buried-source  operator,  can  be  derived 
which  relates  the  seismogram  received  by  a  buried  receiver 
to  that  received  by  a  surface  receiver  lying  on  the  same 
path.  This  operator  depends  only  on  the  layers  lying  between 
the  buried  receiver  and  the  surface  and  could  be  used  to 
relate  the  signals  received  by  sensors  in  a  vertical  array. 

Consider  a  sensor  located  in  a  trivial  boundary  which 
terminates  the  upper  layer  set,  as  shown  in  Figure  7. 


SURFACE 

t*s 

Figure  7.  Buried  receiver 


For  a  wave  X Q ( z )  coming  from  below  and  a  wave  P(z)  returning 
from  above,  the  observed  seismogram  X ( z )  will  be  the  sum  (all 
are  displacement  waves),  i.e. 

X  ( z )  =  Xq(z)  +  P(z) 
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By  (12b) 


p  =  z  M 1 /z )  x 

A  /  _  \  0 


and  from  (7*) 


X 


wkn(l-r. ) 

- —  X 

A  ( z ) 


Combining  we  get 

X  =  A(2)  t  2kA(l/2)  x  (20) 

wkn(l-ri)  s 

This  differs  from  (15)  in  that  the  buri ed-recei ver  operator 
(20)  is  symmetri c.  The  discussion  associated  with  (9)  still 
applies,  that  is,  we  should  not  expect  to  be  able  to  go  from 
the  ghosting  filter  for  our  array  back  to  the  reflection 
coefficients  of  the  layers.  Inspection  of  the  derivation  lead¬ 
ing  to  (11),  however,  shows  that  a  similar  result  does  not 
hold  for  the  symmetric  case,  i.e.  the  symmetric  combination 
is  not  invariant  under  an  inversion  of  the  layers. 

As  with  the  buri ed-source  operator,  we  can  remove  the 
restriction  that  the  sensor  be  located  within  a  homogeneous 
layer  and  allow  it  to  be  located  at  a  non-trivial  interface. 

This  case  is  somewhat  less  interesting  than  the  corresponding 
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extension  for  the  source  operator,  but  it  is  easy  to  do.  The 
result 


X  .  A(2)  +  (1-C)  2kft(l/.)  X 
vrn(l-r.)(l-c)  s 


(21) 


if  the  interface  containing  the  receiver  has  reflection 
coefficient  c. 
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ABSORPTION  AND  THE  FREQUENCY-DOMAIN 
CALCULATION 


So  far  we  have  viewed  our  solutions  in  the  time  domain, 
but  we  could  equally  well  view  them  in  the  frequency  domain. 
The  z-transform  formulation  gives  us  the  vantage  point  of 
Janus  in  being  able  to  see  both  worlds  at  once.  This  is  the 
power  of  the  equal-travel-time  assumption.  While  it  is  not 
strictly  necessary  to  retain  this  assumption  in  the  frequency- 
domain  calculation,  it  is  convenient  for  our  purposes  to  do 
so.  Considering  simplicity,  speed,  and  accuracy,  the  time- 
domain  solution  is  undoubtedly  the  correct  approach  in  the 
absence  of  absorption;  but  since  we  now  wish  to  treat  the 
general  case  of  frequency-  and  depth-dependent  rbsorption, 
a  frequency-domain  approach  will  henceforth  be  necessary. 

(See  Trorey,  1962,  for  an  example  of  what  is  encountered 
in  attempting  a  time-domain  solution  for  this  case).  The 
transition  to  the  frequency  domain  is  easily  made  by  noting 
the  definition  z  =  e“iu)T.  Then  instead  of  representing 
polynomials  in  z,  our  expressions  can  be  considered  to  be 
complex  functions  of  frequency.  Since  we  wish  a  discrete, 
sampled  time  function  with  sampling  interval  t  and  will  be 
using  a  Discrete  Fourier  Transform  (DFT)  to  obtain  it  from 
the  frequency  function,  we  will  take  as  our  phase  interval 
Awt  =  2tt/M ,  where  M  is  the  number  of  time  or  frequency  points 
desired  (see  Rader  and  Gold,  1969,  for  a  discussion  of 
Discrete  Fourier  Transform  theory,  and  Sherwood  and  Trorey, 
1965,  for  a  discussion  of  the  z-transform). 
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Propagator  with  absorption 

Following  Sherwood  and  Trorey,  we  can  obtain  the  layer 
matrix  from  (2)  by  noting  that  z  represents  the  two-way  delay  in 
crossing  the  layer.  With  absorption  present,  the  wave  will 
be  further  modified.  If  we  represent  our  (one-way)  absorption 
law  by  f(z),  then  replacing  z  by  zf  (z)  and  w  by  wf(z)  gives 
the  layer  matrix  with  absorption 


u 

1 

zf2U) 

zf2(z)r 

U 

D 

j  wtf(z) 

r 

1 

D 

J  +  l 


(22) 


We  then  write  the  propagator  (1)  in  the  form 


U 

1 

F'(z) 

G*  (z) 

U 

(23) 

D 

surface 

G(z) 

F(z) 

D 

k 

_  _ 

- 

- 

L  J 

where  F*  and  G'  are  related  to  F  and  G.  (We  take  Claerbout's 
form  for  our  layer  matrix,  rather  than  Sherwood  and  Trorey's. 
In  ours,  the  absorptive  band  f^  lies  above  the  boundary  r^). 
Again,  we  will  be  interested  only  in  the  combinations 
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A ( z )  =  F(z)  -  G * ( z ) 

A'{z)  =  F'(z)  -  G(z) 

defined  by  analogy  with  (4).  We  develop  recursions  for  A  and 
A'  in  the  same  manner  as  in  the  lossless  case.  Multiplying 


the 

propagator  by  an 

additional  layer  matrix. 

F'(z) 

G'(z) 

zfk<z> 

2fk<2>rk 

G(z) 

F(z) 

• 

rk 

1 

- 

_ 

_ 

zfi  F'  +  rkG'  zf l  rkF'  +  G ' 

2fk  G  +  rkF  2fk  rkG  +  F 

from  which  we  obtain 

Ak(z)  =  [F(z)  -  G'(z)]  +  zfk(z)rk  [G(z)  -  F'(z)] 
AJ(z)  =  rk  [G'(z)  -  F ( z ) ]  +  zf^(z)  C F 1  ( z )  -  G ( z ) 3 

or 

Ak^  =  Ak-l^z^  "  zfk^rkAk-l  ^z^ 

Ak(z>  =  ■  rkAk-l<z>  +  zfk(z)  Ak-l(z) 


(24) 


(25) 


corresponding  to  (5).  (To  avoid  conflict*  we  now  use  sub¬ 
scripts  to  indicate  the  iterative  step.)  In  the  absence  of 
absorption,  fk(?)  s  and  i*  1S  easy  t0  show  from  (25)  that 
in  this  case  A‘(z)  =  zkAk(l/z),  which  gives  back  (5).  At  each 
frequency,  (25)  represents  a  set  of  coupled  recursions  over 
layers  in  the  complex  quantitites  A,  A1,  z,  and  f.  Initially, 
A1 (z)  =  A] (z)  =  1 . 

Since  we  require  a  real  time  function  A(t),  A(w)  must  be 
conjugate-symmetric  about  zero  frequency.  This  implies  that  we 
need  only  calculate  A(w)  for  positive  frequencies  (up  to  the 
folding  phase,  -rr ) .  This  also  implies  that  f(z)  must  be 
conjugate-symmetric  (the  "crossing  symmetry"  relations  of 
Futterman,  1962)  in  order  that  (25)  will  lead  to  a  conjugate 
symmetric  function. 

We  will  take  as  our  absorption  law 


fk(w)  =  e’la,Tl/4Qk  +  1<|>k(w)  (2fi) 


where  ( to )  represents  the  dispersive  contribution  to  the  phase, 
and  Q  conforms  to  the  usual  definition  of  the  quality  factor. 
This  "linear  absorption  law"  has  been  experimentally  verified 
In  rocks  (see  Trorey,  1962,  for  references)  and  appears  to  be 
the  best  choice  for  a  calculation  of  this  type.  Neglecting 
dispersion  (discussed  below),  f  =  exp  (- 1 cox | /4Q)  and  we  can 
use  efficient  recursions  over  frequency  to  calculate  the 
exponential  and  the  sines  and  cosines  in  (25). 

We  see  here  the  difficulty  with  the  time-domain  calculation 
of  (25).  For  each  layer  added,  we  must  do  two  complete 
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3 

convolutions.  Thus  this  procedure  is  N  .  For  small  absorption, 
our  filter  will  have  short  effective  length,  so  we  might  still 
get  by.  The  alternative  is  to  use  Trorey's  approach.  In  either 
case,  one  is  also  faced  with  the  necessity  of  obtaining  a  time 
representation  of  (26)  for  each  layer. 


Absorptive  solutions 


Knowing  how  to  calculate  A,  we  can  immediately  transcribe 
all  our  solutions  to  include  absorption  by  the  technique  used 
to  obtain  (22).  We  list  the  results  here  for  reference. 


X  =  wkTT  (l-ri  )fi  (z)/A(z) ;  half-space  (27) 

k  k 

X  =  w  TI[  T  (l-r.)f.(z)  /  Aj(z);  far-surface  (28) 

i  1 

l  k  * 

X  =  wKTT  (l-r.)f.(z)  B*(z)  /  At(z);  bu ri ed-source  (29) 

i  1  1 

X  =  B*(z)  /  wkTT  (l*r. )f^ (z) ;  b.  source  operator  (30) 

X  =  A  ( z )  /  w  TT  (l-r.)f.(z);  b.  receiver  operator  (31) 


To  avoid  conflict  we  have  used  B  to  refer  to  the  distant 
layer  stack,  which  we  have  also  inverted  (stacked  away  from 
observer)  so  as  to  be  consistent.  A  careful  calculation  of  the 
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buri ed-source  operator  for  this  case  shows  (29)  and  (30)  to  be 

the  correct  forms.  Here  A  e  A(z)  -  z  f k ( z )  A'(z),  which  is 

calculated  by  addition  of  the  far  surface  rk  =  1  (which  must 

not  be  included  in  the  transmission  factors  in  28  and  30). 

* 

Evidently  A  is  no  longer  antisymmetric.  Similarly, 

A+  =  A(z)  +  zfk(z)  A' (z). 

Proof  of  minimum-phase 

Claerbout  gives  an  inductive  proof  based  on  (5)  that 
A(z)  must  be  minimum  phase  in  the  lossless  case.  (See  Sherwood 
and  Trorey  for  a  discussion  of  the  terms  "minimum  phase", 
"realizable",  and  "positive-real"  and  their  application  to 
seismic  problems).  By  a  slight  modification  of  his  argument, 
it  is  possible  to  show  directly  from  (25)  that  A(z)  must  in 
fact  ue  minimum  phase  for  any  realizable  absorption  law  f(z). 
We  rewrite  (25)  in  the  form 


Ak(z) 


1  *  zrkfk 


z) 


Ak-l <z> 


(32) 


Ak-i<z> 


zf^(z) 


Ak-l(z) 

Ak-l<z> 


(33) 


Assuming  f(z)  has  no  poles  inside  the  unit  circle  (i.e.  is 
realizable),  then  we  need  only  to  show  that  |  A '  /  A  |  <^1  on  the 
unit  circle  to  complete  Claerbout's  proof.  (It  is  evident  from 
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the  structure  of  (25)  that  A 1 ( z )  must  be  realizable  when  f 
is. ) 


Lemma :  |A'(z)/A(z)|  <  1  when  |z|  =  1 


Proof:  We  prove  (34)  by  induction.  Writing  otk( 
then  dividing  (33)  by  (32)  gives 

zfk  -  rk/“k-i 

“k =  “k-i  ,  ~  - 

1“zrkfk“k-l 

We  must  show  that  |ak|  <_  1  when  |ak.-|l  £  1 » 

I rk  *  2fk  “k-1 I  i  I  1  -  zrkfk“k-ll 

given  that  |z|  =  1,  | f k |  <  1,  and  |rk|  <  1.  We 
by  contradiction,  i.e.,  assume 


rk  -  Zfk“k-1 
1  ‘  zrkfk“k 


.  .2 

ztk“k-l 


'  zrkfk  “k-l 


or 


(34) 

z)  s  AJ(z)/Ak(z), 


1 

. ,  we  require 

(35) 

establish  (35) 
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2  u 

lzfk  ak-l I  +  r2  -  2  Re(rzf2a) 

>  rk I zf k  ak-l I2  +  1  “  2  Re(rzf2a) 


Since  |  z |  3  1 ,  we  have 


ifk  v-iT  +  -k » rkifk  vil‘ + 1 


If  I f k  ak  1  I  =  1  *  we  have  the  desired  contradiction;  if  not, 
we  can  only  have  |f2  ok_1 |  <  1 .  by  the  assumptions.  Collecting 
terms  we  get 


k-1 


2 

]  >  1 


k-1 


2 


Since  the  common  factor  must  be  positive,  we  can  divide  and 
obtain  r2  >  1 ,  a  contradiction.  Since  a-j  =  1  »  the  lemma  is 
provided. 

Having  proved  the  crucial  addition,  we  restate  Claerbout's 
proof.  Assume  Ak_-|  is  minimum  phase.  Assuming  f  is  realizable, 
(32)  must  be  realizable.  Further,  the  second  term  must  have 
magnitude  less  than  or  equal  to  unity  on  the  unit  circle,  by 
the  Lemma.  Hence  its  real  part  must  be  less  than  or  equal  to 
unity  and  therefore  (32)  must  have  positive  real  part  on  the 
unit  circle.  Since  all  Ak  represent  real  time  functions,  (32) 
must  be  real  when  z  is  real.  Therefore  (32)  is  "positive-real" 
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and  hence  minimum  phase.  Thus  Ak  is  minimum  phase  whenever 
Ak_i  is.  Since  A-|  =  1 ,  all  Ak  must  be  minimum  phase. 

The  minimum-phase  condition  on  A(z)  is  necessary  in  order 
that  we  may  obtain  realizable  seismograms.  Inspection  of  (27) 
shows  that  if  our  absorption  law  f(z)  is  also  minimum  phase, 
then  we  actually  obtain  minimum-phase  seismograms,  provided 
we  remove  the  initial  delay.  This  is  true  even  in  the  case 
of  two  free  surfaces  (rk  =  1)  (note,  however,  that  this  may 
introduce  poles  or  zeros  on^  the  unit  circle). 

Pi spersion 

The  key  requirement  in  the  above  discussion  was  the  physical 
realizability  of  f(z).  The  zero-phase  absorption  law  f  = 
exp  ( -  |  ojt  |  /  4Q )  is  evidently  not  realizable,  since  real  symmetric 
frequency  functions  possess  real  symmetric  time  transforms. 

Hence  f(t)  possesses  non-causal  precursor  before  time  t  =  0. 
Futterman  (1962)  has  derived  a  relation  between  the  amplitude 
and  phase  parts  of  (26)  based  on  causality  (Kramers-Kronig 
dispersion  relations).  He  shows  that  the  real  part  of  the  index 
of  refraction  can  be  related  to  the  imaginary  part  by  the 
Hilbert  transform 


Re  An(w)  =  - 


i«"An(n)  dn 
n-w 


(36) 


and  he  obtains  a  simple  form  for  the  phase  in  the  case  of  a 
linear  absorption  law  (26).  He  uses  the  infinite  Hilbert  trans¬ 
form,  but  since  we  are  performing  a  discrete,  frequency-limited 
synthesis,  we  might  expect  the  finite  Hilbert  transform 
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to  work  better.  Such  Is  the  case.  In  (37)  ffl  represents  the 
absorptive  part  of  (26)  and  is  the  dispersive  phase.  (Note 
that  in  this  section  we  use  w  to  mean  wt.  Obviously  the  argu 
ments  in  (37)  must  be  dimensionless.)  It  can  be  shown  (Rader 
and  Gold,  1969)  that  functions  for  which  the  phase  is  the 
Hilbert  transform  of  the  log  magnitude  are  actually  minimum 
phase.  Thus  we  make  our  absorption  law  minimum  phase. 

Putting  (26)  into  (37)  gives 
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(38) 


for  our  phase.  Reducing  (38)  to  the  range  of  positive  fre¬ 
quencies  and  integrating  by  parts  gives  a  term  which  vanishes 
plus  a  term  involving  j  log  sin  u  du,  for  which  no  explicit 
form  exists.  It  can,  however,  be  expressed  in  terms  of 
"Lobatchevsky's  function",  for  which  a  series  expansion  is 
given  by  Gradshteyn  and  Ryzhik  (1965).  Carrying  out  the 
algebra  gives 
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which  is  exactly  the  result  one  would  obtain  by  deriving  the 
frequency-domain  Hilbert  transform  (i.e.  interchanging  sine 
and  cosine  coefficients  and  changing  the  parity)  of  the 
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function  |ft|/4Q.  Thus  we  appear  to  be  stuck  with  (39).  Actually, 
it  turns  out  that  the  most  successful  procedure  is  just  to  do 
the  frequency-domain  Hilbert  transform  using  the  FFT  (Fast 
Fourier  Transform).  Synthesizing  <J>  from  (39)  using  the  inverse 
FFT  was  slightly  less  successful,  supporting  our  original 
supposition  that  when  doing  discrete  synthesis  it  is  best  to 
remain  entirely  within  the  discrete  realm  (observe  that  (39) 
represents  an  infinite  sum  and  the  coefficients  1 / k ^  were 
obtained  by  doing  an  integral).  Futterman's  phase  did  quite 
a  bit  worse.  The  criteria  in  all  these  cases  was  the  observed 
realizability  of  f(t). 

Several  properties  of  the  discrete  dispersive  phase  are 
immediately  evident  from  (39),  namely  4> ( 0 )  =  ( it )  =  0.  Further¬ 

more,  4>  has  an  extremum  at  tt/2  about  which  it  is  symmetric: 

00 

«»(J)  -  -  (-Dk/(2k+l)2  «  -  G/ttQ  (40) 

k=0 

where  G  2  .91596559...  is  Catalan's  constant  and  G/tt  s  .291  5609... 
Using  these  properties,  an  empirical  approximation  to  (39)  was 
obtained , 

<j)(o))  2  -G  log  [1  +  co ( tt-uj )  ]  / ttQ  log  (1  +  tt  /4)  (41) 


While  this  form  was  not  used  in  the  programs,  it  might  be  useful 
in  some  applications.  It  fits  fairly  well,  but  actually  it  is 
not  much  more  trouble  just  to  compute  the  Hilbert  transform. 
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Having  a  procedure  to  compute  the  phase,  we  have  all  the 
elements  needed  to  complete  our  frequency  synthesis.  Using  a 
dispersive  absorption  law  in  (25)  requires  the  explicit  compu¬ 
tation  of  all  the  sines  and  cosines  and  consequently  more 
computing  time.  For  comparison,  routines  were  written  with  and 
without  dispersion  in  A ( z ) .  In  both  cases,  dispersion  was 
included  in  the  straight-through  transmission  filter  (the 
numerator  in  27-29).  This  is  just  given  by  (26)  with  a  dis¬ 
sipation  factor  X -j  1  / Q i ,  which  can  be  large  for  many  layers,  even 
in  the  case  of  high  Q.  Comparison  of  the  semi -di spers i ve  and 
complete  dispersive  calculations  is  given  in  Figure  8.  Running 
time  for  the  dispersive  calculation  was  only  -1.5  times  as  long 
as  for  the  semi -di spers i ve  case,  but  the  semi -di spers i ve  routines 
have  also  been  included  here  because  of  their  simplicity. 

Inversion  methods 

We  have  delayed  until  now  a  discussion  of  the  chief  difficulty 
with  the  frequency-domain  calculation:  aliasing,  as  it  is  called  by 
most  authors.  A  function  that  is  synthesized  in  the  frequency  domain 
and  is  not  time-limited  will  have  a  time  transform  in  which  later 
times  are  folded  over  into  earlier  times,  since  the  computed  spectrum 
is  the  spectrum  of  the  entire  record  and  since  the  DFT  treats  time 
functions  as  though  they  were  periodic  (Rader  and  Gold).  Several 
authors  have  mentioned  the  effect  in  connection  with  this  problem, 
but  evidently  none  has  diagnosed  the  real  difficulty.  The  key,  it 
turns  out,  lies  in  the  method  of  performing  the  deconvolutions 
in  (27-31).  From  time-domain  theory,  we  know  that  the  inverse 
wavelet  A(t)  has  a  finite  length  equal  to  the  number  of  layers 
(in  the  lossless  case),  whereas  the  seismogram  X(t)  has  in 
principle  infinite  length  (and  in  the  case  of  two  free  surfaces, 
i.e.  (28)  and  (29),  does  not  decay  at  all).  This  indicates  that 
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it  would  be  better  to  transform  the  numerator  and  denominator, 
Il.f.(z)  ariu  A(z),  into  the  time  domain  and  do  a  time-domain 
deconvolution  than  to  divide  in  the  frequency  domain  and  then 
transform.  Using  the  first  procedure  we  should  get  zero  aliasing 
error  in  the  lossless  case,  provided  we  take  at  least  as  many 
frequencies,  M,  as  layers,  K.  In  fact,  comparing  the  output 
from  the  time  and  frequency  routines,  LAYERS  and  FLAYERS,  with 
infinite  Q  and  random  reflection  coefficients,  gave  agreement 
to  8-11  significant  digits.  Tests  performed  using  frequency- 
domain  deconvolution,  however,  gave  results  whose  accuracy 
depended  critically  on  M.  The  error  in  this  is  "wraparound  error", 
due  to  the  fact  that  frequency-domain  deconvolution  is  circular, 
whereas  time-domain  deconvolution  is  not.  Thus  the  real  villain 
is  wraparound,  not  al’  "ing. 

In  the  absorptive  case,  A{t)  is  not  strictly  finite  in 
length,  but  develops  a  rapidly  decaying  tail.  Thus  some  aliasing 
error  will  develop  in  this  case,  but  in  general  A(t)  is  still 
niuch  better-behaved  than  X  { t ) ,  and  one  can  make  the  effect 
negligible  by  choosing  M  large  enough.  Since  the  tail  on  A 
decreases  monoton i cal ly* ,  one  can  always  determine  by  inspection 
in  any  case  whether  aliasing  has  been  significant.  This  is 
definitely  not  the  case  with  X.  Values  of  M  about  20%  larger 
than  K  have  been  used  successfully  for  reasonable  Q's,  and  it 
would  probably  never  be  necessary  to  use  values  as  high  as  2K. 


♦Inspection  of  (25)  shows  that  the  highest-order  term  in 
A  ( z )  is  zlcrkn^fi2(z),  so  that  the  tail  behaves  like  the  overall 
transmission  filter  with  twice  the  dissipation  factor.  Thus 
for  any  "buffer"  allowance  M-K,  the  error  will  be  a  function 
of  the  total  dissipation  factor. 
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Another  point  that  should  be  mentioned  here  concerns  the 
way  (22)  was  factored.  If,  Instead,  the  inverse  absorption  were 
left  inside  the  matrix,  then  instead  of  (25)  we  would  have 


Akft>  -  fk'(t)  Ak.,(z)  -  zfk(2)rkAk-l(z) 


Ai<2>  =  -rkfk1(2)  Ak-1(2>  +  2fk<2)Ak-l(2) 


and  in  the  numerators  of  (27-29)  we  would  have  only  a  scalar 
transmission  factor.  This  approach  is  tempting  because  of  its 
simplicity,  but  unfortunately  it  is  numerically  very  unstable 
for  low  Q's.  The  approach  we  have  used,  however,  appears  to 
be  completely  stable  for  the  first  three  cases  (27-29),  even 
for  very  low  Q's.  The  last  two  cases  (30,  31)  are  stable  for 
reasonable  Q's,  but  in  all  likelihood  the  "burying  operators" 
would  be  computed  only  for  rather  few  layers  anyway. 
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EXAMPLES 


Examples  of  calculated  seismograms  are  given  in  Figure  8. 
Nonabsorpti ve ,  absorptive,  and  absorptive-dispersive  seismo¬ 
grams  are  compared  for  each  of  four  different  realizations  of 
random  layer  structures.  Each  structure  consists  of  100  layers 
terminated  at  each  end  by  free  surfaces.  Reflection  coefficients 
were  drawn  from  a  uniform  distribution  on  (-.3,  .3)  and  Q's 
were  drawn  from  a  "log-uniform*1  distribution,  lying  between 
32  and  320,  centered  on  100.  The  source  (a  spike)  was  buried 
10  layers  down  from  the  far  surface.  This  model  is  admittedly 
not  very  realistic  and  is  intended  primarily  to  illustrate 
the  effects  of  absorption.  Two  effects  are  to  be  noted:  a 
general  smoothing  of  the  seismogram,  due  to  the  low-pass 
transmission  filter  in  the  numerator,  and  a  progressive 
simplification  down  the  record,  due  to  the  inclusion  of 
absorption  in  A(z).  Inclusion  of  dispersion  in  A  had  rather 
little  effect,  although  this  naturally  depends  on  the  dissipa¬ 
tion  (39). 
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Figure  8.  Computed  seismograms 

3  '1 


APPLICATIONS 


Of  several  applications  of  the  routines  which  have 
already  been  made,  one  will  be  described  here:  testing  the 
reliability  of  pP  depth-determination  procedures.  Another 
application,  generating  data  to  test  a  homomorphic  filter 
program  designed  to  reduce  convolutional  noise  by  beaming 
in  the  pseudo-time  domain  (complex  cepstrum),  is  being  made 
by  P.  R.  Lintz.  Quasi-real isti c  model  studies  are  also  being 
carried  on  by  R.L.  Sax,  Some  work  has  also  been  done  on 
statistical  behavior  of  synthetic  seismogram  spectra,  which 
will  be  reported  in  the  future. 

Source  depth  determination 

The  numerator  of  the  buried-source  solution  (14)  contains 
the  source  depth  information.  One  might  hope  to  retrieve  this 
information  by  an  analysis  of  the  spectral  zeroes  of  the 
seismogram,  since  the  denominator,  which  is  the  total-path 
filter,  contributes  only  poles.  As  previously  mentioned,  this 
has  been  done  (Cohen,  1970).  The  usual  procedure  is  to  smooth 
out  the  poles  and  enhance  the  visibility  of  the  zeroes  t>y 
averaging  spectra  from  several  stations  for  the  same  event. 

Actually,  the  receiver-end  effect  can  in  practice  be 
factored  out  of  (14).  Assuming  only  weak  reflectors  deep  in 
the  mantle  (disregarding  core  phases),  (5)  says  that  A(z) 
will  be  fairly  short.  The  derivation  leading  to  (15)  applies 
to  any  decomposition  of  the  path,  so  we  may  split  the  path 
in  the  middle.  Then  writing  A  for  the  layers  above  the  source, 
B  for  all  those  on  the  source  end,  and  C  for  the  receiver-end 
layers,  we  can  write  (14) 


X ( z )  -  A*(z)  /  [ B ( z )  C(z)  -  zkT  B ( 1 / z )  C(l/z)] 


I 


If  both  B  and  C  are  effectively  short  compared  with  the  path 
length  ky,  and  if  we  confine  our  attention  to  the  first  part 
of  the  seismogram,  then 


X  ~  A* ( z )  /  B ( z)  C(z) 


Averaging  over  stations  gives 


X  ~  [ A *( z )  /  B( z) ]  •  X-j  Cl/C1  (z)] 


(42) 


assuming  that  A  and  B  are  nearly  the  same  for  each  station. 

The  C.j  will  usually  be  quite  different  and  the  average  will 
tend  to  smooth  out  the  effect  of  these.  Thus  we  are  usually 
justified  in  ignoring  receiver-end  effects  provided  we  average 
over  stations. 

Confining  our  attention  to  the  numerator  A*,  in  the 
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absence  of  reflectors  this  is  just  1-z  for  a  source  buried, 
k  layers  deep.  The  power  spectrum  is  then 
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which  has  nulls  at  integer  multiples  of  fQ  =  1/kt.  Thus 
measuring  the  null  periodicity  gives  the  depth  of  the  source 
in  time  units:  td  =  1 /2f Q .  In  practice,  one  is  usually 
restricted  to  measuring  the  first  null  frequency  because  of 
the  effect  of  the  source  window.  Another  means  of  determining 
the  periodicity  is  to  take  the  spectrum  of  the  spectrum  (the 
cepstrum;  see  Cohen*,  1970). 

Cohen  has  applied  the  method  to  real  data,  with  mixed 
results.  Conversion  of  pP  delay  times  to  depths  requires  the 
application  of  an  "average  velocity"  and  the  resultant  depth 
estimates  will  evidently  be  only  as  accurate  as  the  velocity 
assumed.  Errors  of  this  nature  could  explain  the  discrepancies 
observed  by  Cohen.  Another  possible  source  of  error,  investigated 
here,  is  the  effect  of  reflectors  on  the  null  frequencies.  It. 
turns  out  that  introduction  of  rather  few  layers  with  reason¬ 
able  reflection  coefficients  can  produce  large  shifts  in  the 
nulls.  In  principle,  the  cepstrum  should  provide  a  more  stable 
estimate  of  the  nulls,  but  in  practice  the  cepstrum  can  also 
become  quite  messy. 

The  effect  of  reflectors  can  be  seen  by  writing  the 
source-burying  operator  A  ,  the  numerator  in  (42),  in  the  form 

A* ( z )  =  [1  -  zkA(l/z)/«(r)]  A(z)  (43) 


♦Actually,  Cohen  zeroes  out  the  negative-frequency  part  of  the 
spectrum  before  transforming.  It  can  br,  si, own  that  this  "cepstrum" 
is  just  the  squared  envelope  of  the  autocorrelation.  This  is 
also  evident  from  his  plots.  Our  "cepstrum",  however,  is  just 
the  square  of  the  autocorrelation. 
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Since  |A(l/z)/A(z) |  =  lt  A*(z)  will  still  have  zeroes  on  the  unit 
circle,  although  they  will  no  longer  be  equally  spaced  if 
A ( z )  -/  1.  Evidently  A(z)  itself  will  contribute  no  zeroes  on 
the  unit  circle,  by  Claerbout's  proof  or  by  consideration  of 
(5).  The  number  of  zeroes,  and  hence  the  average  spacing,  will 
therefore  remain  the  same.  Thus  the  cepstrum  should  provide 
a  stable  estimate  where  inspection  of  null  spacing  fails. 

(Large  subsurface  reflectors  would,  however,  put  additional 
zeroes  close  to  the  unit  circle,  which  would  contribute  to 
the  cepstrum  in  the  form  of  peaks  at  shorter  delay  times). 

This  is  shown  in  Figure  9,  in  which  the  ideal  case  is  compared 
with  cases  having  inhomogeneities  above  the  source.  (The  first 
plot  in  this  and  succeeding  figures  is  for  a  homogeneous  medium 
between  the  source  and  surface.  The  other  theee  are  realizations 
of  random  structures,  all  with  the  same  source-depth  delay  time. 
The  random  reflection  coefficients  in  the  latter  cases  are 
drawn  from  a  uniform  distribution  on  [-.3,  .3]  with  zero  mean. 
There  are  15  interfaces  above  the  source,  excluding  the  surface, 
i.e.  this  could  represent  a  delay  time  of  1.6  seconds  with  a 
sampling  rate  of  20  sps.  These  and  most  of  the  following  plots 
are  4-decade  semi-log  plots.)  Although  the  null  spacings  vary 
widely,  the  true  delay  times  (marked  by  arrows)  are  given  in 
each  case  by  the  "break"  in  the  cepstrum.  If  we  consider  the 
denominator  in  (42)  also,  and  allow  layers  below  the  source 
as  well,  then  this  criterion  breaks  down  and  we  see  cepstral 
peaks  at  longer  delay  times  (Figure  10:  reflection  coefficients 
above  the  source  are  the  same;  in  addition,  there  are  15  below 
drawn  from  the  same  population).  Including  a  source  causes 
further  complications  (Figure  11).  (The  source  function  is  an 
empirical  one  of  Cohen's).  Convolution  with  a  source  results 
in  a  multiplication  (or  modulation)  of  the  cosine  ripple  by 
the  source  spectrum.  In  the  cepstrum,  the  result  is  the  source 


-40- 


STRUCTURE  4  1 
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cepstrum  plus  the  echo  cepstrum  plus  the  sum  and  difference 
(beats).  Taking  the  logarithm  of  the  spectrum  before  trans¬ 
forming  makes  the  source  and  echo  additive  and  produces  the 
sum  of  the  two  log  cepstra  upon  transformation.  This  tends  to 
simplify  the  cepstrum  somewhat,  and  appears  to  improve  the 
accuracy  (Figure  12). 

Actually,  the  behavior  of  the  log  cepstra  observed  in 
Figure  12  was  borne  out  by  47  random  realizations  of  the  same 
type,  that  is,  there  was  almost  always  a  sizable  peak  at  the 
correct  delay  time,  although  it  was  not  always  the  biggest 
one.  (In  some  cases,  choosing  the  biggest  peak  would  lead  to 
considerable  errors.)  It  was  expected  that  the  presence  of 
noise  would  tend  to  nullify  any  advantage  possessed  by  the  log 
cepstrum  over  the  linear  cepstrum.  In  a  rough  attempt  to  model 
this  situation,  log  cepstra  were  computed  for  the  same  47 
realizations,  with  the  spectra  clipped  at  1/30  of  the  maximum 
(~  -  15  db).  The  first  three  cases,  with  the  standard,  are 
shown  in  Figure  13.  As  anticipated,  any  superiority  largely 
disappears.  Only  14  cases  of  the  47  yielded  the  correct  answer, 
compared  with  29  correct  answers  in  the  absence  of  noise  (these 
numbers  are  subjective).  In  addition,  there  is  little  tendency 
towards  any  peak  at  all  at  the  correct  delay  time,  much  like 
the  behavior  of  the  linear  cepstra  (which  also  gave  14  correct, 
or  nearly-correcty ,  answers  out  of  47).  The  characters  of  the 
"noisy"  log  cepstra  and  the  linear  cepstra  are  remarkably 
similar. 

Since  our  "cepstra"  are  really  just  the  square  of  the 
autocorrelation,  one  might  suppose  the  autocorrelation  itself 
to  be  less  confusing  by  a  factor  of  two,  since  it  contains 
phase  information.  (Cohen  uses  the  product  of  the  autocorrelation 
with  his  square-envelope  autocorrelation. 
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This  results  in  a  peaky,  smooth  representation,  convenient  for 
visual  analysis,  but  of  course  it  is  no  more  accurate  than  the 
autocorrelation  itself.)  In  Figure  14  we  plot  the  signed  square 
of  the  autocorrelation  for  the  same  cases  so  far  considered. 

(For  convenience  in  plotting,  the  peak  at  zero  delay  is  zeroed  out. 
These  are  linear  plots.)  The  correct  delay  time  should  be 
marked  by  a  large  negati ve  peak.  In  fact,  in  one  realization 
there  are  nj)  large  negative  peaks'.  This  is  of  course  embarrassing, 
and  tends  to  indicate  the  general  unreliability  of  auto- 
correlation/cepstrum  analysis. 

Assuming  then,  that  the  ultimate  arbiter  must  be  the 
position  of  the  first  spectral  null,  since  for  shallow  sources 
it  is  usually  the  only  one  visible,  it  was  decided  to  try  to 
determine  just  how  variable  this  criterion  might  be.  A  hundred 
realizations  of  the  type  discussed  (15  layers  above  the  source 
and  15  below,  reflection  coefficients  between  -.3  and  +.3, 
source  depth  fixed  in  time  units)  were  drawn  and  the  standard 
deviation  of  the  observed  first  null  position  from  the  correct 
answer  determined.  This  turned  out  to  be  about  25%  of  the 
correct  answer.  The  maximum  deviation  observed  is  somewhat 
open  to  debate.  There  are  several  cases  among  the  hundred 
where  it  seems  unlikely  that,  in  the  presence  of  noise  and 
other  factors,  the  first  null  would  be  seen  and  correctly 
identified.  Two  such  cases  are  shown  in  Figure  15.  In  all 
likelihood,  the  second  nulls  would  be  the  ones  picked  here 
(remember  that  these  are  1 ogari thmi c  plots).  Thus  it  seems 
that  this  criterion  could  be  off  by  a  factor  of  at  least  two 
in  some  cases.  Studies  made  with  a  deeper  source  indicate  that 
it  might  be  off  even  more,  the  effect  being  that  a  deep  source 
can  masquerade  as  a  shallow  one.  Three  such  cases  out  of  15 
are  shown  in  Figure  16  (linear  plots).  The  sources  are  4  times 
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Figure  12.  Log  cepstra. 
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Figure  14.  Signed  square  autocorrelation. 
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Figure  15.  Two  interesting  cases 
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as  deep  as  before,  63  interfaces  above  the  source  (drawn 
from  the  same  population  as  before)  and  none  below.  A 
fairly  broad  null  falsely  indicating  shallow  depth  occurs 
in  all  of  these  and,  interestingly,  the  cepstrum  seems  to 
substantiate  the  false  nulls.  The  question  here  is  whether 
the  high-frequency  ripple  in  the  spectrum  would  be  visible 
enough  to  tip  off  the  analyst.  At  least  one  deep  earthquake 
has  been  seen  with  such  a  spectrum  (Cohen,  private  communi¬ 
cation). 

Inclusion  of  absorption  might  be  expected  to  weaken  the 
nulls,  but  for  the  cases  computed  here,  it  actually  had 
rather  little  effect  (Figure  17).  (All  layers  assumed  to 
have  Q  =  50.  The  long  transmission  path  through  the  earth 
was  ignored,  since  it  would  not  effect  the  low  frequency  nulls 
and  since  it  is  more  or  less  included  in  Cohen's  "source  .) 
Non-vertical  incidence,  which  is  not  treated  in  this  paper, 
could  be  expected  to  produce  further  complications. 

In  conclusion,  inhomogeneities  of  the  type  which  may  be 
found  in  sedimentary  regions*  can  cause  large  errors  in  source 
depth  delay  times  determined  from  the  first  spectral  null. 
While  cepstral  analysis  overcomes  part  of  the  difficulty, 
other  effects  can  produce  an  extremely  complicated  cepstrum, 
difficult  to  interpret  and  inaccurate  in  itself.  Thus  this 
method  of  source  depth  determination  would  appear  to  be  of 
limited  utility,  although  the  appearance  of  spectral  nulls 
indicating  shallow  depth  might  be  taken  as  corroborative 
evidence  in  the  presence  of  other  information. 

♦See ,  for  example,  Clark  (1966).  A  velocity  step  of  2.1  to  3.9 
km/sec  gives  a  reflection  coefficient  of  0.3.  The  conclusions 
reached  here  are  not  intended  to  be  rigorous,  but  indicative. 
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APPENDIX 


PROGRAMS 


Three  packages  are  included:  the  lossless,  time-domain 
routine  LAYERS;  the  semi-dispersive,  absorptive,  frequency- 
domain  routine  FLAYERS  (I);  and  the  dispersive  frequency 
routine  FLAYERS  (II).  The  first  package  comprises  the 
routines  LAYERS,  RCTOA ,  and  CONVOLV/POLYDIV.  The  second, 
FLAYERS  (I),  RCAFTOA,  and  CONVOLV/POLYDIV.  The  third, 

FLAYERS  (II),  RCAFTOAD ,  and  CONVOLV/POLYDIV.  The  same 
routine  CONVOLV/POLYDIV  is  used  for  all  three  packages. 

Usage  of  the  routines  is  explained  in  the  introductory 
"comment"  statements.  The  two  frequency  packages  have 
purposely  been  made  interchangeable,  the  only  difference 
in  usage  being  that  a  storage  equivalence  is  allowed  in 
one  that  is  not  allowed  in  the  other.  The  time  routine 
convolves  with  a  source  function.  With  the  frequency 
routines,  the  user  may  do  this  externally,  using  CONVOLV, 
if  he  desires. 

The  subroutines  are  written  in  FORTRAN-63,  a  programming 
language  of  the  CDC  1604-B  computer.  Non-standard  external 
symbols  appearing  in  the  absorptive  routines  are  ERASE  and 
COOL.  Calling  ERASE  (N , X )  zeroes  N  elements  of  array  X. 

Calling  COOL  (LN,  X,  SIGN)  Fast-Fourier  transforms  the 
complex  array  X,  to  frequency  if  SIGN  =  -1.0,  to  time  if 
SIGN  =  +  1.0,  where  LN  =  log2  (number  of  complex  elements  in 
X) (Claerbout  et  al.,  1966).  Overall  execution  time  is  propor¬ 
tional  to  K*M  [or  to  (K  +  KP)«M,  if  IOPT  =  3],  On  the  1604, 
1.25  sec  is  required  for  LAYERS,  12,5  sec  for  FLAYERS  (I), 
and  18.0  sec  for  FLAYERS  (II),  for  K»M  =  104.  Less  time  is 
required  if  some  of  the  reflection  coefficients  are  zero. 
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